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We consider two, apparently similar, models of biological evolution which have been claimed to 
exhibit self-organized critical behaviour. A careful reanalysis of these models, including several new 
analytic results for one of them, suggests that they are qualitatively different. We demonstrate the 
limitations of the mean field description of these systems. We argue that a more precise definition 
of self-organized criticality is desirable and establish several criteria in this connection. 
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O 1 I. INTRODUCTION 

<d : 

The phenomenon of "self-organized criticality" (SOC) has recently emerged as a topic of considerable interest with 
\Q . potential applications ranging from the behaviour of sand piles and the description of the growth of surfaces to generic 
descriptions of biological evolution . It would appear that there is no general agreement on a suitable definition 
of self-organized criticality. (As an example, see the recent discussion in [9|JlO|].) Equally important, there do not 
appear to be universally accepted "black box" tests for its presence or absence based solely on observables. (Formal 
definitions based on microscopic, structural analogies can be useful in model studies. They would seem to have limited 
value in deciding whether any given physical system is or is not self-organized critical.) As a consequence, systems 
with a wide range of characteristics have all been designated as "self-organized critical" . The purpose of the present 
0^ i work is two-fold. We shall point out some striking observable differences between two "self-organized critical" models 
which have a remarkable structural similarity. We will argue that these tests suggest that one of the two models 
should not be regarded as critical. At the least, these differences will suggest that definitions should be refined to 
permit a distinction between such models. Along the way, we shall present a number of new analytic results. 

We shall concentrate on two closely related and simple models introduced by Bak and Sneppen to mimic biological 
evolution. They were introduced in |ll|]l2[ | and further analyzed in Jll| . The models (as we shall use them) involve a 
one-dimensional array of TV sites. Initially, each site is assigned a "barrier" drawn at random on the interval < 6» < 1. 
■ At each update of the model, one identifies the minimum barrier and assigns a new random barrier to this minimum 
site. In addition, we assign new random barriers to a certain number of other sites. (The new barriers are again 
drawn uniformly on the interval [0, 1].) In the first version of the model, we choose the two nearest neighbours of 
' the minimum (the local or nearest neighbour model) In the second version, new barriers are assigned to K — 1 
, additional sites which are chosen at random (the random neighbour model). We adopt periodic boundary conditions 

k> ' ln eac ^ m °del- 

rS \ As we shall indicate, the nearest neighbour model is the richer of the twofl It has a number of fascinating properties. 

The most important of these is that it displays "punctuated equilibrium" |l4],|l5| in which any given (small) segment of 
the sites will experience long periods of inactivity punctuated by brief periods of violent activity. (Hence, its interest as 
a model of evolution.) In equilibrium, the nearest neighbour model has a non-trivial distribution of barriers f p"TJ . One 
finds that virtually all of the barriers are uniformly distributed in the interval 2/3 < hi < 1. In the thermodynamic 
limit, N — > oo, only a vanishing fraction of barriers are in the region < 6,; < 2/3. There is also a non-trivial 
distribution of the minimum barriers selected at each update. The probability of finding a minimum barrier with 
the value b decreases approximately linearly^] to zero as b increases from to 2/3. The fraction of minimum barriers 
having values greater than 2/3 also vanishes in the thermodynamic limit. 



1 The biological interpretation of this model is simple: Each site represents a species. The location of each species represents 
its order in a "food chain". The barrier assigned to each species is a measure of its "survivability". At each update, the least 
survivable species undergoes a mutation and induces mutations in random sites of its nearest neighbours in the food chain. 

2 We believe that this model clearly displays self-organized criticality (in the limit N — > oo). 

3 One can show analytically that this probability distribution is not exactly a linear function even in the limit N — > oo. It is, 
however, very close to linear. In this regard, see Fig. 2 of Ref. Iuj|. 
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One can define an "avalanche" in this model by first establishing a threshold (A < 2/3) and, at each update, 
counting the total number of "active sites" with barriers less than A. An avalanche of temporal extent T is said to 
occur when there are active sites for precisely T consecutive updates of the model. In the limit A — * 2/3, one finds 
(numerically) a power law distribution of avalanches. Specifically, Pa,vai(T) ~ T -108 for large T M|. This power law 
provides the first suggestion of critical behaviour in the model. 

One can also study a variety of correlation functions (following Q|) associated with the space or time movement 
of the minimum barrier. For example, the probability that the minimum barriers at two successive updates will be 
separated by x sites is given by a power law, P x (x) ~ x~ 3,1 for large x [p"l| . The most natural temporal correlation 
to study is the probability that, if a given site was the minimum at some time, it will also be the minimum at time 



t later. For large t, this correlation is also given by a power law, Pt(t) ~ t~ 0A2 |17|. Power-law correlations in both 
space and time would seem to be the most compelling indicators of the presence of a critical phenomenon. We argue 
below that power laws in temporal correlation functions alone are probably not a sufficient criterion for criticality and 
that subsidiary conditions on the exponents are needed. 

To make this point, consider the simpler random neighbour variant of the model jlj] in the thermodynamic limit. 
This model also has a non-trivial distribution of barriers. In this case, virtually all barriers are uniformly distributed 
over the interval 1/K < hi < 1. The fraction of barriers in the region < hi < 1/K vanishes. There is again 
a non-trivial distribution of the minimum barriers selected at each update. The probability of finding a minimum 
barrier with the value b is constant over the region < b < \[K. Only a vanishing fraction of the minimum barriers 
have values greater than 1/K. It can be shown analytically]^ that avalanches again follow a power law distribution 
(provided that one chooses the threshold to be A = 1/K) with P avaj i(T) <~ T~ 3 / 2 . 

The spatial correlations of the random neighbour model are completely different from those found in the original 
local neighbour model. In the random neigbour model, there is a significant probability that the minimum barrier will 
be at the same site at two successive updates. As could be anticipated from the fact that there is no notion of spatial 
ordering in the model to begin with, there are no other spatial correlations in this model. Spatial correlations have 
been degraded from power laws to a delta function. However, the temporal correlation of the minimum (as defined 
above) is still given by a power law, Pt(t) ~ i -3 / 2 for large t (see Sec. ^v|) . We will argue in Sec. that the existence 
of this power law alone is not a good indicator of criticality. The value of the exponent is crucial. We will argue that 
the normalizability of Pt(t) provides a useful indication that the system is not critical. By this criterion, we would 
conclude that the random neighbour variant does not describe a critical system. The local neighbour variant does. 

The original motivation for the random neighbour model was precisely that it was easy to solve. In particular, it 
allows for a straightforward and relatively successful mean field description of the barrier distribution [jl2|. However, 
the simplification process results in the degradation of those power-law correlations which provide the strictest test of 
criticality. In fact, this mean field description is not particularly successful at describing the most interesting features 
of the model which are related to the distribution of barriers. The most interesting "observable quantity" to be 
extracted from the barrier distribution is the avalanche probability distribution. For this we need the distribution of 
"active" sites with barriers below 1/K which (by definition) contribute to avalanches. These barriers comprise only 
a vanishing fraction of of all barriers in the thermodynamic limit. Further, all avalanches must end. (Otherwise, 
there would be little merit in talking about avalanche distributions.) Thus, the fluctuations in the number of active 
barriers must be comparable to its average value. This contradicts the usual assumption of mean field theory that 
fluctuations of quantities are small compared to their average value. Under such circumstances, the nature of a mean 
field description renders it of limited utility. It is good at describing the features of the critical state of the system. It 
is not suitable for providing a description of the interesting fluctuations about this critical state. On the other hand, 
the number of sites with barriers larger than 1/K is O(N). In this case, fluctuations about the mean field results 
are genuinely small (i.e., of 0(\/~N)), and their mean field description can be (and is) successful. Unfortunately, such 
sites are largely inert and of little physical interest. 

The organization of this paper is as follows. In Sec. || we review the mean field approximation to the random 
neighbour model [|l2| and an exact analysis of this model in terms of the number of active sites JD| . We show how 
to obtain the distribution of barrier heights exactly and contrast the results with the approximate values obtained 



with mean field theory. In Sec. Ill we shall indulge in a brief digression in order to demonstrate that the power 
laws in "avalanche" distributions are not necessarily indicative of critical behaviour. We shall make this point by 
considering a one-dimensional random walk and an equivalent one-dimensional chain of non-interacting spins. (We 
take the position that these systems should not be regarded as critical. Others may disagree.) In Sec. IV we consider 



the role of spatial and temporal correlations in the random and local neighbour models. Here, as noted above, we 



4 The case K = 2 has been discussed in [h3| , and the general proof is given here in Appendix 
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will encounter a significant difference between the random and local neighbour models. A number of conclusions will 
be drawn in Sec. |V[ along with suggestions for possible further tests of criticality. Two appendices are provided for a 
detailed account of the derivation of various analytic results. 



II. MEAN FIELD THEORY AND BEYOND 

In this section we wish to summarize the description of the random neighbour model using mean field theory [ fL2[ . 
Since it is our intention to assess the validity of mean field theory, we will also introduce a simple and tractable 
reformulation of the model ]l3[ | which leads to exact results without the need for simulations. 

In the mean field theory approach to the random neighbour model, we replace the joint probability distribution 
p(bi, . . . , 6jv) of the barrier heights by p(b) N , where p{b) is the average equilibrium distribution of barrier heights at 
one site and N is the number of sites in the model. It is easiest to work with 4>(b) defined as 

0(6) = f db' p(b') . (1) 

Jb 

Evidently, p(b) follows immediately by differentiation. Further, the probability distribution C(b) for the value of the 
minimum barrier is also determined by 4>(b), 

By demanding equilibrium one can derive the following simple equation for 4>{b), 

(N -K)(j> N + N(K-1)4> + K(N -1)(6-1) = . (3) 

Here, we have allowed K, the total number of barriers changed at each update, to assume any value K > 1. 

Now, let us rewrite the random neighbour model by concentrating on the number of active sites, A, at any given 
update |l3| . An active site will be any site which has a barrier height less than some threshold value < A < 1 . The 
number of active sites can change at any update of the model. The probability for a given change will depend on 
the current number of active sites. A detailed analysis of this active site description is performed in Appendix ^ (for 
finite N) and in Appendix |b| (in the thermodynamic limit N — > oo). 

To obtain a sense of the nature of this model, it is useful to consider the average number of active sites as a function 
of the threshold A. To do this, consi der t he average change in A in one update if A > 1. The calculation is done in 
Appendix |A|, and we obtain [see Eq. (A5)] 



AiM^-D- ^f;" ■ <4) 

First, consider the case A < 1/K. The value of AA is less than for all values of A and, more importantly, for 
all values of N. As we approach the thermodynamic limit N — > oo, we see that AA will approach some limiting 
negative number with additional negative corrections of 0(1/N). It is clear that the average value of A A will be 
positive when A is precisely zero. Thus, for A < 1/K, the average number of active sites has some value, A(X), which 
evidently depends on A but is independent of N. The density of active sites, defined as p = A/N, is precisely in the 
thermodynamic limit for all A < 1/K. 

For A > 1/K the situation is quite different. Now AA is positive for all values of A less than some Nx c (X) and 
negative for all A greater than Nx c (X). The value of x c (X) is readily seen to be 

K — 1 N K — 1 y ' 

In this case, the average number of active sites will grow with N. The density of active sites will be x c (X). In other 
words, in the thermodynamic limit, 

_ f for A < 1/K 

P ~ [ (KX- 1)/(K- 1) for A > 1/K ' W 

The density of active sites grows linearly from to 1 as A grows from 1/K to 1. At the point A = 1/K, the density is 
continuous but its derivative is discontinuous. This result is superficially suggestive of a second-order transition with 
a critical point at A = 1/K. 
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Let us now determine the average number of active sites for A < 1/K in the thermodynamic limit J^] In Appendix 0, 
we have derived the generating function for the probability of finding precisely A active sites, f2(z) = Xm>o ^ > (^)- s 



Z' 



One finds 



n(*) = (i-/'(i)) (1 „ f {z) . (7) 
/(«) - z 



where 



/(z) = (l-A + Az) K ■ (8) 

From this, the explicit values of P(A) can be extracted by means of contour integration p3[ . The average value of 
the number of active sites is 



A = £ = = ^A 1 -^ 1 ^ 2 . (9) 



As noted, this result applies for A < 1/K. The root-mean-square deviation from A can also be calculated from Sl(z). 
We find that 



5 a = fn"(i) + - (fi'(i)) 2 ] 1/2 



VKX 
1 - KX 



1 - |(isT + 1)A + + l)(5Jf + 4)A 2 - ^-K(K + 1)(K + 5)A 3 
2 D 12 



"I 1/2 

(10) 



The important point to note is that 5 A is larger than A for all A < 1/K and equals A at the limit. 

The result of Eq. (||) is interesting. It allows us to see short-comings of the mean field analysis of this model for 
A < 1/K. In the mean field approach, 

Ami (A) = N[l - <f>(\)} . (11) 
It is easy to solve the mean field approximation in the limit of large N and for A < 1/K. We find 

(1 - K/N)(l - A mi /N) N + (K- 1)(1 - A m{ /N) + K(l - 1/N)(X - 1) = (12) 

which leads to 

A mf = -ln(l-ifA) (13) 

in the thermodynamic limit. While Eqs. (|l3| ) and (JsJ) agree in the limit A — * 0, they differ everywhere else. These 
differences are particularly striking as A approaches the critical value of 1/K. 

To obtain a quantitative sense for the differences between mean field theory and the active site model, we compute 
the equilibrium distributions of barriers, p(X), and of the minimum, C(A), in each approach for finite N and compare 
the results with simulations. In the mean field approximation, we have to solve Eq. (^) numerically for 0(A). Using 
Eqs. (|l|)-(||), we can derive explicit expressions for p(X) and C(A) in terms of 0(A). In the active site version, p(X) is 
obtained from 

1 dA , . 

pW = nTx- ^ 

Note that this definition agrees with that of mean field theory through Eq. (|l|). Since A — Xm>o P{A)A, we obtain 

A>0 

where P'(A) denotes the derivative of P(A) with respect to A. The probability distribution of the minimum is simply 



5 This quantity follows immediately from Eq. (^) when A > 1/K. 
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C(A) = -P'(0) 



(16) 



in the active site description. In Appendix |X] we describe how to obtain P{A) and P'(A) for finite N. Numerical results 
for p(X) and C(A) obtained from mean field theory, the active site analysis, and simulations for K = 3 and N = 100 
are displayed in Table Q. The barrier and minimum distributions for the active site and mean field calculations are 
again different. As expected, the (exact) active site results are in agreement with the results of simulations. The 
barrier distribution in the region A < 1/K is, by nature, small. Thus, absolute differences between mean field and 
active site results are also small. However, this region contains the active sites which carry the interesting physics. The 
fractional difference between mean field and active site results is large, and we regard this disagreement as significant. 

Of course, the mean field results are sometimes completely reliable. In the thermodynamic limit and for A > 1/K, 
we can solve the mean field equations by noting that <p is now genuinely less than 1. This justifies dropping the term 
<f> N to obtain exactly the result of Eq. (|5|) obtained with the active site calculation. This success is not surprising. In 
this domain, A grows linearly with N while 8 A grows like v^V. 

It is easy to understand these results. The success of a mean field calculation of any quantity requires that the ratio 
of the rms deviation of this quantity to its average value vanishes in the thermodynamic limit. When the deviation 
of a quantity is of the same order as its average value, the mean field approximation is likely to fail. Indeed, the fact 
that 5 A > A for all A < 1/K is a clear indicator that mean field approximations are not likely to be reliable in this 
region. Since one anticipates large fluctuations to be the hallmark of any critical system, there is virtually an a priori 
contradiction in expecting that mean field theory will provide a complete description of any critical system. 



III. POWER LAWS WITHOUT CRITIC ALITY 



Some authors have used the existence of power laws in the distribution of avalanches as a function of their duration 
as an indicator of criticality in similar systems. We have already mentioned in Sec. | that both the local and the 
random neighbour model show such power laws with exponents —1.08 (for the local version with A = 2/3) and —3/2 
(for the random version with A = 1/K). One notes the existence of avalanches "of all size scales" (meaning no 
exponentials) and declares criticality. It is thus useful to consider a system which is non-critical by common consent 
and which nevertheless has such power laws — a random walk in one dimension. 

In this case, we define the start of an "avalanche" as being when the walker is at the point x = 0. Now consider 
a (P-type) walk for which the walk ends as soon as x is again 0. We are interested in P 2 „ which is the probability 
that a walk will return to x — (for the first and only time) at step 2n. To determine this probability, it is useful 
to consider a Q-type walk of length In for which x = at the end points without regard for whether x — is also 
reached at intermediate steps. This walk occurs with probability Q2n- Evidently, Qi n is trivial to determine from 
combinatorics: 

Q2n = ^( 2 ") ■ (17) 



2 2n \ n 

Any Q-type walk is either a pure P-type walk or can be decomposed into some sequence of P-type walks of varying 
lengths provided only that the sum of their lengths is 2n. Using this observation, we deduce that the generating 
function for the numbers Qm is given by 



1 + Q^Z 2n = oo 

n=l 2^n=l^nZ 



(18) 



Summing over all walks has led to great simplification. The various terms contributing to the right hand side at a 
given order merely describe the various ways in which a Q-walk can be built from P-walks. 

Fortunately, the left side of Eq. ( |l8| ) can be summed analytically given the form of Eq. (|l7|). Thus, 

oo 

1 - P ^ z2n = Vl-^ 2 ■ (19) 

n=l 

This allows us to write 

_(2n-3)H [2 1 
P2n ~ (2n)H V n (2^7* (20) 

where the last expression is the large n limit . 
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It is no accident that an identical power is found for the avalanche distribution of the random neighbour model. 
The results of the last section (and of the next) indicate that this distribution function also comes from a random 
walk of a very similar character. A f-step avalanche starts when the number of active sites changes from zero to 
non-zero and ends when the number of active sites is again at step t. The rules for changing the number of active 
sites, given in Appendix pi are simply a set of random walk r ules. The probability for an avalanche of length t (with 



A = 1/ K ) in the random neighbour model is found from Eq. ( |B14[) to be 



aval 



(t) = J^^- 3/2 + <?(t- 5/2 ) • (2D 



This exponent is identical to the random walk exponent be cau se the problems are essentially identical. This is only 



true for X — l/K since then the average change AA of Eq. (A5) equals zero (in the thermodynamic limit) so that we 
are dealing with an unbiased random walk. 

To underscore the fact that this particular power law is not a good indicator of criticality, consider a one-dimensional 
chain of spins ±1 which do not interact.^] As usual, we adopt periodic boundary conditions. For simplicity, we consider 
only the thermodynamic limit as the number of sites goes to infinity. The only question is how to define an avalanche. 
We choose to consider the distribution of the sizes of domains of zero magnetization. In other words, start at any 
site and count the number of sites until the total magnetization is zero. The first zero encountered defines the size of 
the domain. The ensemble average distribution of domain sizes, P(t), is precisely given by the random walk problem 
solved above. Thus, P(t) ~ i~ 3 / 2 . Clearly, this power law does not indicate that the non-interacting spin lattice is 
critical. 

Since there is at least one example of a manifestly non-critical system which exhibits a power law in an avalanche 
distribution, the presence of such power laws cannot be regarded as a sufficient indicator of critical behaviour. 



IV. THE ROLE OF CORRELATIONS 



So far, the behaviours of the local and random neighbour models are qualitatively identical. They each have a 
non-trivial distribution of barrier heights and a corresponding distribution of the heights of minimum barriers. Each 
has a power law avalanche distribution. It is our a priori expectation that the local neighbour model is a genuine 
self-organized critical model. The fact that the random neighbour variant is essentially identical to a random walk 
and, hence, to a non-interacting one-dimensional spin chain strongly suggests that this model is not critical. We would 
thus like to find a test which distinguishes between these models. It would seem natural to consider the spatial and 
temporal correlation between localized avalanches. Since a sound definition of local avalanches is somewhat subtle, 
we consider instead the spatial and temporal correlations between the minimum barrier (following Q). This has the 
virtue that it is very much like the kind of question one would ask to determine criticality in condensed matter systems 
or in statistical mechanics. It is also a sensible first step. For the local neighbour model, all sites in a single, local 
avalanche must have been activated during the avalanche. For every site in the local avalanche, the minimum barrier 
must have passed through the site itself or through one of its nearest neighbours. Hence, the path of the minimum 
provides the "skeleton" upon which the full avalanche will be built. Another virtue of working with the minimum 
barrier is that there is neither need nor room for the somewhat artificial introduction of a threshold or any special 
definition of an active site. 

In the nearest neighbour model, results for the spatial and the temporal correlation functions can currently only 
be obtained by simulations. The spatial correlation function P x {x) is defined as the probability that the minimum 
barriers at two successive updates are separated by x sites. There are two temporal correlation functions of interest. 
The first-return probability, S{t), is defined as the probability that, if a given site is the minimum at time £o ? it will 
again be the minimum for the first time at time to + t. The all-return probability, Pt(t), is the probability that this 
site will also be the minimum at to + t regardless of what happens at intermediate times. In all cases, we observe 
power law behaviour for large arguments. Specifically, for large x we have @ 

• (22) 



This system is not critical even when there is an interaction. Everything which follows can be generalized to include an 
arbitrary nearest-neighbour interaction and a magnetic field of arbitrary strength using text book techniques. 
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For large t the first-return probability is found to be 



The all-return probability for large t is |jj 



S(t)~^ • (23) 



m-^M ■ ( 24 ) 



Numerical results for the spatial correlation function and the temporal correlation function (first- and all-returns) 
were obtained by simulations and are shown in Fig. [j] and Fig. ||, respectively. The spatial correlation function is 
familiar from Bak and Sneppen pdj | who describe the motion of the minimum as a random walk in which the distance 
moved by the minimum goes like?/ 3 . For any number of sites N, the all-returns temporal correlation function follows 
the power law of Eq. ( pi[ ) until it hits the accidental coincidence rate of 1/N. It does this when 

At^0.05N 2A . (25) 

The interpretation of this result is simple. The power law correlations do not extend over times longer than the 
longest (local) avalanche in the system. The spatial extent of this longest avalanche will be on the order of N. The 
above equation tells us roughly the relation between the largest (spatial) avalanche and its (temporal) duration. Note 
that the exponent 1/2.4 = 0.42 is not equal to the value of 1/3 which one might crudely have expected from Bak and 
Sneppen. While it remains to be seen, we might expect that the duration of an avalanche of spatial size n will be of 
order n 2A . 

Let us now consider the random neighbour model. Here, the correlations can be obtained analytically. First, consider 
the spatial correlation between the minimum at successive updates. There is some probability that the minimum will 
stay at the same site. In the thermodynamic limit, this probability is given by pi = 1 — (K — 1) \n[K/(K — 1)], see 
Eq. ( B30[ ). Otherwise, the distribution of next minima is random (i.e., spatially uniform). Results for this correlation 



function have also been obtained by simulation and are shown in Fig. || As expected, there is complete agreement 
between the analytic results and the simulation. 

Analytical results for the temporal correlation function are more difficult to obtain. An exact calculation of the 
first-return probabilities, S(t), is given in Appendix [b] in the thermodynamic limit. Similar calculations for finite N 
are possible but tedious. The leading term in S(t) in the thermodynamic limit is 



s w = ir*^'~ 3/2+ ° (t ~ B/2) • (26) 

Additional terms in the asymptotic expansion of S(t) are given in Appendix |b| Convergence of S(t) for finite N 
to the values given here in the thermodynamic limit is remarkably slow. Thus, Fig. |J shows a comparison of these 
analytic results with the first-return probabilities arising from simulations for the case N = 10 5 . While the agreement 
is extremely good, it is not perfect. Since the arguments in Appendix |b] are exact, several comments are in order. 
First, the small discrepancies apparent for small t are not due to finite N effects. Rather, they are due to our use of 
(three) leading terms an the asymptotic expansion of S(t). In spite of this truncation, agreement is excellent. The 
discrepancies at large t are finite N effects. 

Analytic calculation of the all-return probability, P t (t) , is more challenging. The naive expectation is that all-return 



probabilities can be determined using an equation similar to Eq. (18), i.c 



In fact, this equation is not exact since it involves ensemble averages at an inappropriately early stage in the calculation 
(see Appendix |B]for details). It is far from trivial to solve for Pt(t) exactly due to reasons which will be discussed in 
Appendix^. There, we derive an improved equation for P*(t). However, an explicit analytic solution of this equation 
is, at best, tedious. Therefore, we have elected to perform an approximate calculation using Eq. (|27|). We can use 
the results of simulations to test the approximate validity of Eq. (^?j) in relating all-returns to first-returns. It is 
surprisingly accurate for both the local and rando m ne ighbour models and results in errors which are, at worst, at 
the 1% level. The result for large t is given in Eq. (|B4l|) in the thermodynamic limit. We find 
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For completeness, we have also performed simulations for this quantity. Again, convergence with N to the thermody- 
namic limit is extremely slow. Thus, the results shown in Fig. ^ are for a simulation with 10 5 sites. The values of Pt(t) 
assume the expected i~ 3 / 2 form quite quickly. This form persists for some time. Since the number of sites is finite, 
there is a finite probability (independent, of time) that a given site will "accidentally" be chosen to be updated. For 
sufficiently large t, this accidental rate of 1/N will dominate Pt{t). This fact is clearly reflected in the form of Fig. |] 
for very large times. The discrepancies for small t are larger than those observed for the first-return probabilities. 
This has nothing to do with deviations from the thermodynamic limit or, worse, inadequacies in the use of Eq. (p^). 
Rather, it is related to the fact that the asymptotic expansion of Pt{t) converges slowly. This discrepancy can be 
removed by noting that, as a consequence of Eq. P*(l) = 5(1). (This result is evidently trivial and exact.) 

The next few terms can be improved analogously by expanding Eq. ( [27| ) rather than using the asymptotic expansion 
Eq. (^8|). Except as noted in the asymptotic region of large t, the discrepancies between our calculated values of Pt(t) 
and simulations are at the level expected from the approximate validity of Eq. ([27]). 

We are now in a position to assess the main qualitative differences between the two versions of the model. The local 
version displays power law correlations in both space and time. We view this as a compelling indicator of criticality. 
The random variant of the model has no spatial correlations apart from a finite same-site correlation. It would seem 
unfair to place too much emphasis on the absence of spatial correlations in the random neighbour model since this was 
clearly built in by design. The random neighbour model does show power law behaviour in the temporal correlation 
function, and it is legitimate to discuss the differences in this correlation function for the two versions of the model. 

We believe that there are two reasons why the power laws of the random neighbour model are not indicators 
of cr itical ity. The first argument is of a qualitative nature. The arguments presented in Appendix M leading to 



Eq. (B37) rely solely on the analogy to a random walk. They require nothing more than the ensemble averages of 



various one-body operators. We believe that a convincing demonstration of criticality must rely, either direc tly o r 



indirectly, on the demonstration of the existence of genuine correlations. Since the arguments leading to Eq. (B41) 
do not require any information about two- or many-body operators, we do not believe this power law (or any other 
power law obtained in a similar fashion) should be regarded as a convincing indicator of criticality. 

The second argument is more quantitative. While Eq. (|27l ) is not generally valid for either the random or local 
models, we show in Appendix [b] that it is true for the random neighbour model in the special case z = 1. This allows 
us to relate the sum over all-return probabilities to the sum over first-return probabilities as 

1 + |* (B) = '-sLsO) • (29> 

It is clear from the definition of the first-return probability that the sum over all S(n) must be finite and cannot 
exceed 1 (J2n>i ^( n ) = £ < 1). By contrast, the sum over all-return probabilities need not be finite since any given 
site can be visited any number of times over an infinite time interval. If S < 1 , Eq. ( p9| ) indicates that the sum over 
the Pt(n) will be finite. When this happens, it means that a given site will be active only a finite number of times 
even though we study the system for an infinite time interval. It is not reasonable to regard as critical any system 
in which individual sites "die". On the other hand, if S is exactly equal to 1, the sum over the Pt(n) diverges. In 
this case, sites will never die. Thus, we believe that a normalized sum over the S(n) (S = 1) and, equivalently, a 
divergent sum over the Pt(n) should be necessary requirements for criticality. (The equivalence of these criteria is 
general and does not require the validity of Eq. (|27|). If, upon reactivation, there is a finite probability that the site 
will die, the sum over all-returns must be finite. In practice, it is usually easier to detect a divergence than it is to 
extract a specific, finite normalization.) Obviously, the condition for convergence or divergence of the sum over the 
Pt{n) can be reformulated in terms of the exponent a in the asymptotic form Pt(n) ~ n a . If a < — 1, the model is 
not critical. 

The random neighbour model does not satisfy these requirements. We see from Eq. ( B35| ) with z = 1 that 



Sn>i S(n) = ]-/K < 1. Thus, at any update when a given site becomes active, there is a probability 1 — 1/K that 
it will never become the active site again. The reason for this is clear from the known distribution of all barriers and 
of minimum barriers: The distribution of all barriers indicates that O(N) barriers have heights greater than 1/K. 
No barrier with height greater than 1/K will ever be the minimum barrier (in the thermodynamic limit). The new 
value drawn for the barrier at the active site has a probability of 1 — 1/K of being assigned a height greater than 
1/K. When this happens, the site dies and will never become the minimum again. As a consequence, the sum over 
Ptin) must converge. In fact, Pt(n) ~ n~ 3 / 2 for large n. Hence, the random model does not satisfy this criterion 
for criticality. These considerations only apply in the thermodynamic limit of infinitely many sites. This provides no 
limitations since it is only in this limit that the question of criticality arises. 

The situation is dramatically different in the local neighbour model. Simulations reveal that Pt(n) ~ n~ 0A2 for 
large n. Thus, the sum over the Pt{n) must diverge, and the sum over first-return probabilities must be identically 
1 in the thermodynamic limit. The local neighbour model passes this test of criticality. On the naive basis of the 
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one-body distributions of all barriers and minimum barriers, this result is surprising. A newly selected site will 
draw a new barrier of height greater than 2/3 with a probability of 1/3 and might be expected to die. Of course, 
this site can be re-activated through the presence of spatial correlations. The value of this barrier can be reset to 
a value less than 2/3 provided only that its nearest neighbour becomes the minimum. The mechanism for site re- 
activation provides the connection with our previous point: One-body arguments cannot provide this re- activation, 
and dynamical correlations are required. 

Power law spatial correlations are a good indicator of criticality. Power laws in temporal correlators are only 
useful indicators of criticality when they meet the additional criteria noted above. Given this and the fact that 
these additional conditions on temporal correlators are normally enforced through the agency of spatial correlations, 
it might seem pointless to consider temporal correlations at all. This is not the case. It is easy to devise sensible 
models for which the extraction of information regarding spatial correlation is either difficult or impossible. For 
example, consider a "quenched" version of the current model in which the additional sites that are to be changed 
when a given site becomes the minimum are chosen according to some complicated (but fixed) scheme according 
to which we guarantee only that each site has the same number of "neighbours" . (In particular, the designation of 
"neighbour" need not be reversible.) In this case, the definition of a spatial separation between sites may be subtle 
and may be ambiguous. However, the temporal correlation function for the minimum can be studied as above even for 
such systems. Temporal correlators can permit decisions regarding criticality even in cases where spatial correlation 
functions are unavailable. 



V. DISCUSSION AND CONCLUSIONS 



We have considered a local neighbour model of evolution and a simplified random neighbour variant. The random 
model was originally introduced to permit an approximate analytic description of the model using mean field theory. 
We have demonstrated that it is potentially dangerous to modify critical models in order to facilitate mean field 
calculations. Mean field theory assumes that fluctuations are small while criticality insists that they are large. In 
practice, mean field theory does not provide an adequate description of the distribution of active sites even in the 
random neighbour model. Such disagreements are not merely quantitative. Consider, for example, the distribution 
of all barriers as A = 1/K — SX approaches the threshold value of \[K. For any finite value of SX it is always possible 
to consider a value of N sufficiently large that the results of Eqs. (|J) and (^) (obtained in the thermodynamic limit) 
apply. These results indicate that p m f (A) ~ 1/5 A. The exact results indicate that p(X) ~ 1/5X 2 . The results of Table | 
indicate the presence of similar qualitative errors in the mean field description of the distribution of minimum barriers. 
Fortunately, as we have indicated here, most of the interesting aspects of the random neighbour model permit exact 
analytic solution. 

Both local and random neighbour models have non-trivial distributions for all barriers and minimum barriers. Each 
possesses a power-law avalanche distribution. However, the T~ 3 / 2 distribution encountered in the random neighbour 
model is precisely analogous to that of a one-dimensional random walk. It does not strike us as a reliable indicator 
of criticality. The discussion in Sec. Ill suggests that it is never safe to draw conclusions regarding criticality on the 
basis of avalanche distributions. 

A more detailed consideration of spatial and temporal correlations between minimum barrier locations reveals a 
clear distinction between the two versions of the model. The local neighbour model displays power law correlations 
in space (with exponent —3.1) and time (with exponent —0.42 for all-returns). The random neighbour model lacks 
the necessary power-law correlations in space (by design) and shows power law correlations in time (with exponent 
—3/2 for all-returns). 

It seems reasonable to extend the definition of self-organized criticality in order to distinguish between these models. 
Thus, we think it useful to require simultaneous power law correlations in space and time. (Definitions of criticality 
which focus more closely on correlations are closer in spirit to those adopted in statistical and condensed matter 
physics. This would appear to be a virtue.) Moreover, we require that the temporal correlation function be non- 
normalizable (i.e., it must not fall off faster than 1/t). This is equivalent to the requirement that sites must never 
"die" in the thermodynamic limit. The local neighbour model passes these two tests while the random neighbour 
model fails. This also suggests that any simplified version of a critical model must, at the very least, retain non-trivial 
correlations in its solution. If all observables in a model can be calculated from one-body operators, it is likely that 
those features of the model responsible for criticality have been approximated away. This appears to be the case for 
the random neighbour model considered here. 

Power laws are ubiquitous in self-organized critical models. They appear in (i) the duration of "global" avalanches, 
(ii) the two-body spatial correlation between minimum barrier sites, and (iii) the two-body temporal correlation 
between minimum barrier sites. There is also a power law which relates the total number of active sites in a "global" 
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avalanche to its duration, T.[] It is probably safe to insist that all four of these quantities must be power law if a model 
is to be considered self-organized critical. Power laws in all of these four quantities are also a necessary condition for 
a space-time plot of the model to be self-similar. It is tempting to speculate that the best definition of self-organized 
criticality may be the self-similarity of its space-time plot. Of course, self-similarity is a stronger statement than the 
mere presence of power laws. This question merits a closer look. 
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APPENDIX A: ACTIVE SITE ANALYSIS FOR FINITE N 

Here, we derive some analytical results for the active site description of the random neighbour model for finite N. 
At each update the value of the barrier at the minimum site and at if — 1 other sites, selected at random, are replaced 
by new barriers which are drawn at random on the interval [0, 1]. 

In general, we start from a configuration in which there are precisely A active sites with barrier heights less than 
some threshold value A. We then consider how the number of active sites can change following one update of the 
system. These update rules result in a description of the system as a random walk in A with probabilities for specific 
changes in A which depend on the current value of A. At each update, the net number of active sites can be reduced 
or increased. The first step in each update is to remove the K barriers which are to be changed. This can reduce the 
number of active sites. If A = 0, no active sites can be removed. If A > 0, we remove the minimum barrier (which 
must be an active site) and r — 1 additional active sites (where < r — 1 < K — 1). This reduces the number of active 



sites by r. The associated probabilities, d r ^Ai for reducing the number of active sites by r are given as 1 1 3 1 



A-l \ f N-A 
r-1 ) [ K-r , 

drA = ~ / \ ■ (Al) 

N — 1 
K - 1 

One assigns K new random barriers which can create a number of new active sites, t, ranging from to K . The 
probability for creating t active sites, Ct, is given by 

c* = ( f ) A'(l - A)*"* - (A2) 

The Ct are just the coefficients of z f in the expansion of f(z) = (1 — A + Xz) K . (We everywhere use the convention that 
= if to < or to > The fact that the location of the K — 1 additional sites is chosen at random implies 

to ' 

that there are no spatial correlations in the system (other than a finite same-site correlation). The coefficients d r ^A 
and Ct provide complete information about the system. In principle, knowledge of their values (for all A) is sufficient 
to calculate all observables in this model. 

Now, let P{A) be the equilibrium probability of having A active sites and let Q(z) — Xm>o P{A)z a be the 
corresponding generating function. The condition for equilibrium is obtained by requiring that the P (A) should not 
change when the system is updated: 

K K 

Q(z) = J2 P{A)z A = P(0)f(z) + J2 P{A)z A ]T d r . A z- T J2 ctz* 



A>0 A>0 r=l t=0 



A>Q 



(A3) 



7 The total number of active sites in an avalanche is roughly the area of the avalanche in a space-time plot. 
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where gA{z) — Y^=i d r .AZ~ r for A > 0. Unfortunately, g depends on A, and it is not possible to deriv e a simple 
analytic expression for Q(z) except in the thermodynamic limit (see Appendix |bJ). Nevertheless, Eq. (A3) leads to a 
set of N homogeneous equations for the P(A) which are soluble by standard numerical techniques for any value of N 
when supplemented by the condition that the sum over all P(A) equals 1. 

It is useful to note that a much faster numerical approach is available in the special but important case when 
A < 1/K. In this case, only a finite number of the P(A) contribute (independent of TV), and the sums in Eq. (A3) 
can be truncated safely. These equatio ns can be solved by a simple iteration starting from, e.g., P(0) = 1. 

Before discussing ways to solve Eq. ( |A3| ) exactly, it is useful to compute an exact result needed in Sec. |l|. If A > 0, 
the probability for a change A A in the number of active sites in one update of the model, Pa A, is 



PAA = ^ d r ^ACAA+r 



(A4) 



r=l 

Assuming that A > 0, the average change in A in one update is 

K-l 



AA= PAAAA=(fg A y(l) = K\ + g' A (l) 

(A-1)(K-1) 



AA=-K 

(KX - 1) - 



N - 1 



(A5) 



Let us now return to the solution of Eq. ([A3]). As noted, the most straightforward approach is to rewrite it as a 
matrix equation for the P(A) which can then be solved numerically. In order to compute the equilibrium distribution 
of barriers and of the minimum we also need a matrix equation for the P'(A) (defined as dP(A)/d\). Using the 
normalization condition ^2 A>0 ^(^) = 1' we obtain for A > 1 



P(A) = 

This leads to the matrix equation 



N 



1 - P(A') 



K 



CA 



d r ,A-t+rC t P(A -t + r) 



r=l i=0 



N 



J2 Maa<P(A') = Y(A) , 



(A6) 



(A7) 



where 



and 



K 



Maa' = 5aa> 



CA 



^ GA-A'+rdr,/ 



r=l 



Y(A) = c A 



(A8) 



(A9) 



To obtain p(A) and C(A), we differentiate Eq. (A6) with respect to A. Note that the only additional A-dependence 
lies in the c^'s. Using Xm>o ^'(^) = 0, we obtain for A > 1 



K K 



P'(A) 



Here, c' t is simply 



-ca Y P '( A ') + C A P (°) + d 7%A -t+r[c' t P(A -t + r) + c t P'(A - t - 

A' = l r=l t=0 



/ 



dc t 



{t- KX)X t - 1 (l- X) 



K-t-l 



We again obtain a matrix equation with the same matrix, £a'=i M aa >P'(A') = Z(A), where now 

K N 



Z{A) = c' A P(0)+J2 E c' A _ A , +r d r>A P{A') 



(A10) 



(All) 



(A12) 



r=l A' 
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The results quoted in Table | were obtained by solving these two matrix equations numerically for K = 3 and N = 100 
and applying Eqs. (|l5| ) and (|l^). 

While this procedure is straightforward, it is not elegant. We observe that one can also use Eq. ( |A3| ) to derive a 
(K — l)-order differential equation for 



Q(z) 



n(z) - o(o) 



The equilibrium distributions of barriers and of the minimum are then obtained from 

p(A) = —— O'(l) 
Py 1 NdX V ' 

and 

Note that the normalization 0(1) = 1 yields O(0) = 1/[1 + 0(1)] 



(A13) 



(A14) 



(A15) 



Using the definition of Eq. ( A13) and the explicit form of gA (z), Eq. (A3) can be rewritten as 

-l K 



zfl(z) 



1 



/(*) 



= 1 



N 
K 



r=l 



A>0 



N-A\ P{A)z A ~ 1 
K-r ) O(0) 



Replacing A — 1 by the operator zd z we obtain 



zh(z) + 1 

m 



1 



N - 1 
K- 1 



-1 K 



zd z 
r-l 



N - 1 - zd z 
K-r 



O(z) 



Here, a binomial coefficient involving an operator T is defined as 

T{T - 1) • • • (T — m ■ 



For example, in the special case of K = 2 one obtains 

zO(z) + 1 



/(*) 



= l + 0(z) + 0'(z) 



1 - z 
N - 1 



(A16) 



(A17) 



(A18) 



(A19) 



Eq. (A17) represents an ordinary linear differential equation for O(z) of order K — 1 . In particular, it allows for a 
systematic expansion in powers of 1/ {N — 1) through an iterative solution of Eq. (A.17). Such an iterative procedure 
is completely algebraic in nature and does not re quire t he sol ution of a differential equation. The quantities O(z), 
p(X), and C(A) can then be determined by Eqs. (A13)-(A15). While the matrix equation above provides explicit 
expressions for general K, it is cumbersome for large N and A > 1/AT where the numerical effort required to invert the 
N x N matrix is appreciable. The method just discussed is more suitable for the computation of finite N corrections 
for any K. 



APPENDIX B: ANALYTICAL RESULTS IN THE THERMODYNAMIC LIMIT 



It is possi ble to obtain a great many more analytic results in the thermodynamic limit (N — > oo). In this case, the 
d r ^A of Eq. (Al) are strictly zero for all r except for di : A = 1- Thus, gA{z) = z^ 1 for all A > 0. This simplification 
permits the removal of gA(z) from the sum in Eq. (|A3|), and we obtain 



O(z) = O(0)/(z) + (0(z) - O(0))z-7W 
Together with the normalization 0(1) = 1 this yields 

,(!-*)/(*) 



Q(z) = [!-/'(!)] - 



m - z 



(Bl) 



(B2) 



with f(z) given after Eq. (A2). 
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1. Avalanche Distributions 



In order to find the probability, q(t), of an avalanche of length t, defined with respect to the threshold A, we must 
look in detail at the random walks defined by z^ 1 f(z). Let \{ z ) = St>o ^W 2 * denote the generating function. 
Further, define q n (t) as the probability for having zero active sites for the first time at the t-th time step, given 
the fact that at time zero there were n active sites. Denote the generating function for q n (t) by Xn( z )- (Since an 
avalanche starts and ends with zero active sites, x( z ) = Xo( z )-) The probability distribution of the number of active 
sites after the first time step is the same whether we start with zero or with one active site at time zero. It follows 
that x( z ) = Xi{ z )- To compute x„, it is crucial that the random walk defined by z _1 /(z) can decrease the number 
of active sites by at most one. Therefore, if we start with n active sites, there will be a unique first time when there 
are exactly n — 1 active sites, a unique first time when there are exactly n — 2 active sites, and so on. The probability 
distribution of these time intervals is precisely given by xi, and therefore 

Xn(z)= X (z) n ■ (B3) 

We are now in the position to derive an equation for %. Suppose that, at time step zero, there is one active site. 
Then, the probability for having r active sites after the first time step equals the coefficient of z r in f(z). If there are 
r active sites at time zero, then the probability distribution of the first occurrence of zero active sites is determined 
by Xr( z ) = xi z Y ■ Combining these facts we see that the probability distribution for the first occurrence of zero 
active sites, having started with zero active sites at time zero, is zf{x{z))- The factor of z is due to the fact that we 
take one "extra" time step at the beginning. On the other hand, this quantity is precisely the avalanche probability 
distribution. Thus, x satisfies 

X = z f(x) ■ (B4) 
In order to find a power series for x-i we need several auxiliary results. First, for any integers A and Y > 0, 

y 



E 



=i 



Al 

Y - 1 I \ I 



(-1)'=0 . (B5) 



The proof is as follows: the first binomial is the coefficient of z l in z^ Y ~ l ^/ A {l — z x l A )~ Y . The second binomial times 
(— 1)' is the coefficient of z Y ~ l in (1 — z) Y . Thus, the total sum is the coefficient of z Y in the product of these two 
expressions. The product is 

(l + z 1 /.4 + z 2/A + ... + z ( J 4-i)/^ y z (y-i )M f (B6) 



but this does not contain z Y so that the sum in Eq. (B5) vanishes. This implies that, for all A 



S 1 -^''-'!^^ 31 (B7) 



by expanding the powers of (1 — A) and rewriting the resulting product of binomials in the form of Eq. (B5). 
Let us now introduce the quantity y defined by 

y^zA(l-A)*- 1 (B8) 

and the function u(y) as the solution of 

ufoMl-ufo)]*- 1 ^ (B9) 
which has the property that u(0) = 0. We claim that 



u(l-X) 

* - W^) ■ (B10) 



Indeed, 
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fix) = (1 - A + \ X ) K 



(1- 


X) K 


(1- 


u) K 


(1 - 


- X) K 


2/(1- 


- u)/u 


u{l 


-A) 


z\(l 


-u) 


X 




z 





Now, we are in a position to prove that the power series for \ is 



X^^l-A)^- 1 )^- 1 



2>0 



Kl 



I J (K- 1)1 + 1 



The proof reads as follows: 



X(z) = '£z l (l-\)( K -V l + 1 \ 1 - 1 



l>0 



1 



Kl 



Kl 

I J (K -1)1 + 1 
1 



I / (K —1)1 + 1 
Kl 



(l-A)A- 1 -„)(*-i 
v ' z>o 



1 



Z / (if — 1)Z + 1 



A7 



i / (if - 1)Z + 1 



_ u(l - A) 
_ \{l-u) ' 

Using this expression, we find for large t that 
1-A / 



q(t) 



K 



A W27r(if-1) 



(1 - X) K - l XK K (K - l)-^- 1 ) <- 3 / 2 + 0(t- 5 / 2 ). 



(Bll) 



(B12) 



(B13) 



(B14) 



The quantity in square brackets is 1 for the "critical" value of A = 1/if. For this value of A, we recover the typical 
random walk power law with an exponent of —3/2. 



2. Temporal Correlations 

The final quantity of interest is the temporal correlation function. Define S(b, A; At) to be the probability that, 
given the fact that at some time a certain site is the minimum with barrier b, it will after At time steps again be the 
minimum, for the first time, with value A. Let A; z) = X^At>o ^(^' At)z At denote the corresponding generating 
function. To evaluate this quantity we note that the value of the barrier at this particular site must already be A after 
the first time step^] After the first time step, a certain number of barriers will have a value less than A. We have to 
wait until they have all disappeared before the site in question is again the minimum. We have already computed 
these probabilities, they are just the coefficients of \r where there are r active sites below the barrier, A. It remains to 
compute the probability P± (b, A; r) that, given the fact that the value of the minimum barrier was b, the new barrier 
value at this site is A, and that the number of active sites with value less than A is r (after the update). Then, we 
have 



8 The possibility of an "accidental" change of the barrier at this particular site (through the other (K — 1) randomly selected 
sites) has probability zero in the thermodynamic limit and can be neglected. 
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(B15) 



r>0 



Let Po{b, A;r) denote the probability that, given the fact that at time zero the minimum is b, there are r sites with 
value less then A. We do not need the explicit form of Pq. It is sufficient to note that, when we integrate over b, we 
will obtain the probability for having r sites with barriers smaller than A: 



dbC(b)P (b,X;r) = P(r) 



(B16) 



From this, we can compute Pi . 

Since the minimum is assigned the value A, the number of sites with value smaller than A cannot have changed by 
more than (K — 1). In addition, we must distinguish between the cases r = and r > 0. (For r = 0, b must be larger 
than A. For r > 0, b must be smaller than A.) This leads to 



Pi(M;r) 



K-l 
r 



E 

t>0 



(l-A)^ r A r F„(M;0) 

A (l-A) /f - 1 -'-+ t - 1 A''- t + 1 Po(M;i) ■ 



K-l 

r-t + 



This yields 



f dbC(b)e(b,\;z) = zy2 



K- 1 
r 



(l-A)*- 1_r A r P(0) 



E 

t>0 



z[l-X + X X (z)] 



K-l 



n(o) 



K— 1— r+t-1 \r— 1+1 

fi(x(X))-n(o) 



p(t) 



= z[l-A + A X (zr- 1 [l-/'(l)] 



1 - X(*) 



/(x(*))-x(*) 



We can further simplify this expression by using the form of x given by Eq. (BIO). We find 

6(A,z) = f db C(b) 9(5, A; z) = z(l — A + Ax)^ 1 [1 - /'(l)] ) ~ X 
Jo f(X)-X 



l-A + A X JVA; /( X )-x 
z(l-tfA)(l-u)x 1-X 



(1-A) zx/z-x 
z(l-KX)(l-u) 

(l-z)(l-A) 1 * J 
z(l- KA)(1 -it) f X-u 



(l-z)(l-A) V^( 1 -'«) 
z{l-KX) it 



(l-z)(l-A) 

To obtain a power series for 0, we need one for u. We claim that 

«(y) = f* {Kb ~ 2) [ 

fe((if-l)6-l)!H 

We have already demonstrated that 



A 



7?/ 



1 - u 



E 



2ft 



1 



I J {K- 1)1 + 1 



y 



(B17) 



(B18) 



(B19) 



(B20) 



(B21) 
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Now, 



E 



(Kb-2)\ b _ K -l r / Kb\ y 



-,y 



E 



UK - 1)6 - lV.br K 

b=l u ' ' b>0 

K - 1 



b Kb -I 



K 2 



-I) 



k I dy y 1 * [ K b ) y ' 

J b>a v J 



where J dy stands for the operation y r — ► y r+1 /(r + 1). From Eq. (B21) we derive that 



E 

b>0 



Kb 

b 



y" = 



(K-l)y— + l 
dy 



1 - u 



(B22) 



(B23) 



By differentiating u(l — u) = y with respect to y we find 



u'(y) = 



u(y)(l - u{y)) 
y(l-Ku(y)) 



Using this result, we find that 



and 



Armed with this we obtain 



E 

b>0 



Kb 

b 



V = 



Ku 
1 - Ku 



d 

dy 



y K u =y 



, i (K-l)u 



K(l-Ku) 



g (Kb-2)\ nib _K-l 



6=1 



{{K- 1)6-1)16!" 



K 2 



y k I d yy 1 ^E( & ) 



K - 1 i 



K 



y' K I dy y 1 " 



h>0 



1-Ku 



= y K 

= u . 



k I dy i( y ^ Ku ) 



This proves the claimed power series for u. 
We have now obtained the result 



6(A;z) 



z{l-K\) { 1- (Kb -2)1 

(i-^i-^y-x^dK-Db-iy.bi t zX{1 - x) 
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The next step is to integrate over A. If A > 1/K, the site will never again become the minimum. Thus, we need only 
integrate from to 1/K. The integrals involved are easy, since 



1 d r 



(1 - KX)X b - 1 (l - A)^- 1 ) 6 - 1 = ~ A b (l - A)^ 1 '" 



bd\ 



(B29) 



The final result is 



1/K 



6(» = / dX 9(A; z) = n - (if - 1) In 
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(B30) 
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Unfortunately, it is not easy to extract the behaviour of the coefficients for large powers of z from this expression due 
to the factor 1/(1 — z). This behaviour can be obtained mo re ea sily from a power series constructed about z = 1 



instead of z = 0. The desired series can be derived from Eq. (B30) by writing it as 



dz 

where u(z) is the solution of 



d ( (l-z)Q(z)\ ^ u(z) 



which satisfies u(0) = 0. The function u has a series expansion 



uil-uf-^z ^ 1 . (B32) 



= J? (l + £%(l-^' /2 ] , (B33) 



j>0 



and from Eq. ( |B32| ) one can easily determine a few of the coefficients Oj. That this is a series in \/l — z rather than 
(1 — z) is due to the fact that the function f(u) — u(l — u) K ~ 1 has a maximum at u = 1/K. Given the series 
Eq. ( p333| ), ®{z) is recovered from 

B{z) = ^-f\^ , (B34) 

where the proper integration limits follow from the requirement that 0(z) should have no singularity as z — ► 1. The 
first few terms in the expansion of O(z) are 



z j i 2 /2(^-1)^ 1/2 , 5JT-4, 



47K 2 - 71^ + 26 / 2(^-1) _ 3/2 

90^(^-1) V X 1 j 
2(134if 3 - 291X 2 + 204X - 46) 2 
+ 405K 2 (if - 1) (1 ~ Z) 

6409AT 4 - 17954JC 3 + 18261 K 2 - 7964^ + 1252 / 2(K-1) 5/2 
15120X 2 (Jf - l) 2 V K ( ~ Z > 

+0((l-zf)} . (B35) 

Note that 9(1) = 1/K as expected. Thus, there is a probability of (K — 1)/K that a site "dies" upon updating and 
will never become active again. 

For K = 2, there exists a simple closed expression for 8(z), 



e(z) = {y/T=z-hi(l + y/T=z)\ . (B36) 

1 — z 



The series expansion, Eq. (B35), for 9(z) around z = 1 determines the asymptotic behaviour of S(At) for large At, 



and we find once more a power law with the (leading) random walk exponent —3/2, 



?(.\t\-J 2{K - l) f 1 rArt-3/3 I 14A ' 2 - 2i ^- 13 rA ^-5/2 

^ (A ' } - V^^\3^ (A<) + mK 2 (K-i) {At) 

169if 4 + 142if 3 - 543X 2 - 80K + 313 7/2 
+ 2016K 3 (K- l) 2 ( ^ 

+G((Ai)- 9/2 )} . (B37) 

These are the "first-return" probabilities of Sec. It is far more difficult to compute the desired temporal correlation 
function P t (At), i.e., the "all-return" probabilities. The naive expectation would be to generalize Eq. ([l8|), leading 
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to Eq. (p7j). This is not correct, however. The underlying assumption is that the ensemble of configurations at the 
moment of a first return (after a fixed number of steps) is equal to the ensemble of all configurations. One can show 
numerically that this assumption is not satisfied. This implies that, in order to perform an exact calculation of the 
all-return probabilities, one needs the simultaneous probability distribution of all barriers below a given threshold at 
all times. Fortunately, substitution of z — 1 in Eq. (^|) corresponds to averaging over all ensembles of configurations 
at the moment of first returns. This means that the previous objections do not apply in this situation and that 



Eq. ( |29| ) is, indeed, exact. This is important for the discussion in Sec. IV. 

One can derive an improved equation for the all-return probabilities by taking into account some of the correlations 
that have been neglected in Eq. (|27]). Specifically, we can account for the different probability distribution of the 
minimum in the two ensembles mentioned previously (at the moment of a first return and averaged over all config- 
urations). The improved generating function for the all-return probabilities, &(z) = X)a*>o Pt(At)z At , can thus be 
written as 

j-l n „i/K 

$( z ) = y2 dAiC(Ai)TT/ d\ l+1 Q(Xi,X i+1 ;z) . (B38) 

n>l Jo i=l Jo 

Note that this generating function can also be obtained from 

,1 ,1/K 

$(z) = / daC{a) / dp ${a,P\z) , (B39) 



o 



where $(a, P; z) is given by the integral equation 



rl/K 

$(a,P;z) = Q(a,P;z) + I dX 8(a, A; z)$(A, P; z) . (B40) 



Eq. (B40) does not appear to admit explicit analytic solution. Thus, we have elected to use the approximate rela- 



tionship, Eq. fl29|), to obtain the asymptotic behavi our o f the all-return probabilities. While it is not simple to invert 



1 — O(z) in closed form, we can again employ Eq. (B3J) to extract the asymptotic behaviour of Pt(At) for large At. 
After some algebraic effort we obtain 



2(K-1) f K /A N , /9 A2K 2 + 24K + 1 , A . 

(At) 



1521K 4 + 802K 3 - M17K 2 + 36K + 17 



-7/2 



18144^(^-1)4 

+G((Ai)- 9/2 )} (B41) 

with the now-familiar leading exponent —3/2. This approximate solution for the all-return probabilities does exceed- 
ingly well in describing the data as can be seen from Fig. 0. 
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FIG. 1. Simulation results for the spatial correlation function of the minimum as defined in Sec. |iv| for the nearest neighbour 
model (iV s itcs = 500, 10 s updates). 

FIG. 2. Simulation results for the temporal correlation function of the minimum as defined in Sec. [Fy| for the nearest 
neighbour model (iVsites = 500, 10 s updates). Crosses indicate first-returns, diamonds all-returns. 

FIG. 3. Simulation results for the spatial correlation function of the minimum for the random neighbour model with K = 3 
(Nsitcs = 500, 10 5 updates). 

FIG. 4. Simulation results for the temporal correlation function of the minimum for the random neighbour model with K — 3 
( A^ite s = 100, 000, 2 x 10 7 upda tes). Crosses indicate first-returns, diamonds all-returns. Also shown are the predictions of 



Eqs. (B37) (solid line) and (B41) (dashed line) 
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TABLE I. Distribution of barrier heights, p(A), and of the minimum, C(A), for N = 100 sites and K = 3 as a function of 
A. Results are shown for the mean field approximation, the (exact) active site model, and the simulation of 10 8 updates. The 
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